#!/usr/bin/env python
import numpy as np
from astropy.io import fits
import sys
import matplotlib.pyplot as plt
from scipy.stats import sigmaclip
from testlib.imgconfig import load_imgconfig
from astropy.table import Table

dark_f = sys.argv[1]
mask_f = sys.argv[2]
chip_type = sys.argv[3]
outd = sys.argv[4]
conf = load_imgconfig(chip_type)
dark = fits.getdata(dark_f)
mask = fits.getdata(mask_f,dtype='>i2')
darks,ratio08s,ratio16s = [],[],[]
for ch in range(conf.nchan):
    x1,x2,y1,y2 = conf.mos_x0[ch],conf.mos_x0[ch]+conf.nx[ch],conf.mos_y0[ch],conf.mos_y0[ch]+conf.ny[ch]
    dark_ch,mask_ch = dark[y1:y2,x1:x2],mask[y1:y2,x1:x2]
    dark_clipped,l,h = sigmaclip(dark_ch.flatten(),5,5)
    dark_mean = np.mean(dark_clipped)
    flag08,flag16 = np.zeros(mask_ch.shape),np.zeros(mask_ch.shape)
    flag08[np.where((dark_ch-dark_mean)>0.08)] = 1
    flag16[np.where((dark_ch-dark_mean)>0.16)] = 1
    flag08[np.where((mask_ch>32))] = 0
    flag16[np.where((mask_ch>32))] = 0
    onum08 = np.sum(flag08 == 1)
    onum16 = np.sum(flag16 == 1)
    print(onum08,onum16)
    darks.append(dark_mean)
    ratio08s.append(onum08/conf.ny[ch]/conf.nx[ch])
    ratio16s.append(onum16/conf.ny[ch]/conf.nx[ch])
print(np.arange(conf.nchan))
tab = Table()
tab['chan'] = np.arange(conf.nchan)+1
tab['dark_current'] = darks
tab['ratio08'] = ratio08s
tab['ratio16'] = ratio16s
tab.write(outd+"/dark_chan.tab",format='ipac',overwrite=True)
#fits.writeto("dark_outlier.fits",data=flags,overwrite=True)
